%--------------------------------------------------------------------------
% computes the eigenvalues of the linearized problem using a 'frozen time'
% approach
%
% assumes the basic state is in the t = O(1 / delta) regime
%--------------------------------------------------------------------------


function [ev, kk, M] = comp_eigs(kk, p, h0)

if nargin == 1
    p = params;
end

if nargin < 3
    h0 = 1;
end

% basic state at z = h_bar
v_0 = (1 - p.beta) * (1 - 1 / h0);
c_b = p.beta + v_0 - 1/2 * p.delta * (1 - p.beta) * (p.beta + v_0) + p.delta * p.beta * (1 - p.beta) / 6;

N = p.N;
z = linspace(0, h0, N);
h = z(2) - z(1);

M = spalloc(N, N, 4 * N);


for j = 1:length(kk)
    
    k = kk(j);
        
    if (kk(j) < 0.01)
        w = @(z) p.Ma * z.^2 .* (k^2 / 4 * (1 - z) + p.delta / p.Ma / 2 * (3 - z));
    else
        A = w_coeffs(1, h0, kk(j), p);
        w = @(z) (A(1) * z + A(2)) .* cosh(k*z) + (A(3) * z + A(4)) .* sinh(k*z);
    end
    c1bz = @(z) -p.delta * (1 - p.beta) * (p.beta + v_0) * z / h0^2;
    
    for i = 1:N
        
        if(i == 1)
            M(i,i) = -2 / h^2 - k^2;
            M(i,i+1) = 2 / h^2;
            M(i,N) = -w(z(i)) * c1bz(z(i));
        elseif (i == N - 1)
            M(i,i) = -2 / h^2 - k^2;
            M(i,i-1) = 1 / h^2;
            M(i,i+1) = 1 / h^2 - w(z(i)) * c1bz(z(i));
        elseif (i == N)
            M(i,i) = -2 * (1 + h * p.delta * (1 - 2 * c_b)) / h^2 - k^2 - w(z(i)) * c1bz(z(i));
            M(i,i-1) = 2 / h^2;
        else
            M(i,i) = -2 / h^2 - k^2;
            M(i,i-1) = 1 / h^2;
            M(i,i+1) = 1 / h^2;
            M(i,N) = -w(z(i)) * c1bz(z(i));
        end
    end
    

    tmp = sort(eigs(M, 5, 'SM'), 'descend');
    ev(j) = tmp(1);
end

if nargin == 1

    
    % Lambda = O(1) all k
%     for j = 1:length(kk)
%         if (kk(j) < sqrt(p.delta))
%             w = @(z) p.Ma * z.^2 .* (kk(j)^2 / 4 * (1 - z) + p.delta / p.Ma / 2 * (3 - z));
%         else
%             A = w_coeffs(1, 1, kk(j), p);
%             w = @(z) (A(1) * z + A(2)) .* cosh(kk(j)*z) + (A(3) * z + A(4)) .* sinh(kk(j)*z);
%         end
%         sigma(j) = -kk(j)^2 + (p.Lambda * quad(@(z) w(z) .* z, 0, 1) - (1 - 2*p.beta)) * p.delta;
%     end
%     loglog(kk, -ev, 'k', kk, -sigma, 'r*', ...
%         kk, (1 - 2*p.beta) * p.delta * ones(size(kk)), 'b--', 'linewidth', 1);
%     xlabel('$k$', 'interpreter','latex','fontsize',12);
%     ylabel('$-\sigma$','interpreter','latex','fontsize',12);
%     l = legend('numerical', 'asymptotic', '$(1 - 2\beta)\delta$','location','best');
%     set(l, 'interpreter','latex','fontsize',11);
    
    % Lambda = O(1 / delta)
%     sigma = -kk.^2 + (p.Lambda / 80 * kk.^2 - (1 - 2*p.beta)) * p.delta + 11 / 40 * p.Lambda / p.Ma * p.delta^2;
%     loglog(kk, -ev, 'k', kk, -sigma, 'r*', ...
%         kk, (1 - 2*p.beta) * p.delta * ones(size(kk)), 'b--', 'linewidth', 1);
%     xlabel('$k$', 'interpreter','latex','fontsize',12);
%     ylabel('$-\alpha(0)$','interpreter','latex','fontsize',12);
%     l = legend('numerical', 'asymptotic', '$(1 - 2\beta)\delta$','location','best');
%     set(l, 'interpreter','latex','fontsize',11);
   
    
    % Lambda = O(1 / delta^2)
%     semilogx(kk, ev, 'k', 'linewidth', 1);
%     xlabel('$k$', 'interpreter','latex','fontsize',12);
%     ylabel('$\alpha(0)$','interpreter','latex','fontsize',12);
%     ylim([-0.1 1/2]);
%     xlim([1e-2, 3]);
%     f = @(x) interp1(kk, ev, x);
%     k1 =  fzero(@(x) f(x), 1e-1);
%     k2 = fzero(@(x) f(x), 2);
%     hold on;
%     plot([1e-2, 3], [0, 0], 'k--', k1, 0, 'r*', k2, 0, 'r*');
end